subroutine readvec(lth, vec)
implicit none
integer lth
real*8 vec(lth)
!f2py intent(in) :: lth 
!f2py intent(out) :: vec(lth) 
open(22,file='vec_cfour.data', form='formatted')
read(22,*)vec(1:lth)
close(22)
return
end

subroutine save_eri_df(naux, p, q, iin)
implicit none 
integer naux, p, q
real*8 iin(:,:,:)
!f2py intent(in):: naux, p, q
!f2py intent(in):: iin(naux,p,q)

integer*8 naux8, p8, q8
integer i, j, k
naux8=naux
!write(6,*)'naux8', naux8
!write(6,*)'naux8', naux
p8=p
q8=q
open(22, file='eri_df_aa', form='unformatted')
write(22) naux8, p8, q8
write(22) (((iin(i,j,k),i=1,naux),j=1,p),k=1,q)
!write(22) (iin(i),i=1,p*q*naux)
close(22)

return
end

subroutine putlstso(nirrep, nfirst, nblock, leth, ispin, &
                    ilist, filename, zlist)
implicit none
integer nfirst, nblock, ispin, ilist, leth
integer nirrep
character*8 filename
real*8  zlist(nblock*leth)

!f2py intent(in) :: nfirst, nblock, ispin, ilist
!f2py intent(in) :: nirrep
!f2py intent(in) :: filename
!f2py intent(in) :: zlist(nblock*leth)
integer*8 irrepx
integer i, j, k
integer*8 popp(8), popq(8), popr(8), pops(8)
integer*8 ioppq(8), ioprs(8), irrep0(16)
integer*8 moionuml(10,20), moionumr(10,20), rectmp(10,20)
character*10 filename2
real*8, allocatable :: tmp(:)
integer iprcln, iprcwd, nwrdf, nwrdf2
integer*8 nwords, nstart, nend, recstart, recend, irecord
integer*8 istart, iout, listrec

iprcwd = 4096
iprcln = 4096*8

filename2(1:8) = filename
filename2(9:10) =  '.c'
OPEN(UNIT=413, FILE=FILENAME2, FORM='UNFORMATTED')
rewind 413
read(413) IRREPX
read(413) (popp(i), i=1,8)
read(413) (popq(i), i=1,8)
read(413) (popr(i), i=1,8)
read(413) (pops(i), i=1,8)
read(413) (ioppq(i), i = 1,8)
read(413) (ioprs(i), i = 1,8)
read(413) (irrep0(i), i = 1,8)
read(413) (irrep0(i), i = 9,16)
read(413) ((moionuml(i,j),i=1,10),j=1,20)
read(413) ((moionumr(i,j),i=1,10),j=1,20)
read(413) ((rectmp(i,j),i=1,10),j=1,20)
close(413)
OPEN(UNIT=412, FILE=FILENAME, FORM='UNFORMATTED', access='direct', recl=iprcln)

listrec = rectmp(ispin, ilist)
nwords = nblock*moionuml(ispin, ilist)
if(nwords<=0) then
   return
endif
nstart = (nfirst-1)*moionuml(ispin, ilist)
nend   = (nfirst+nblock-1)*moionuml(ispin, ilist)

recstart = (nstart-1)/iprcwd
nwrdf = nstart - iprcwd*recstart
recstart = recstart + listrec + 1

recend = (nend - 1)/iprcwd
nwrdf2 = nend - iprcwd*recend
recend = recend + listrec + 1

allocate(tmp(iprcwd))
if(recend>recstart) then
   read(412, rec=recstart) tmp
   do i = 1, iprcwd-nwrdf
      tmp(nwrdf+i)=zlist(i) 
   enddo
   write(412, rec=recstart) tmp

   do iout = 1, recend - recstart - 1
      irecord = recstart + iout
      istart = (iout-1)*iprcwd + iprcwd - nwrdf + 1
      write(412, rec=irecord) zlist(istart:istart+iprcwd-1)
   enddo

   read(412, rec=recend) tmp
   istart = nwords - nwrdf2
   do i = 1, nwrdf2
      tmp(i) = zlist(istart + i)  
   enddo
   write(412, rec=recend) tmp
elseif(recend==recstart) then
   read(412, rec=recstart) tmp
   do i = 1, nwrdf2 - nwrdf
      tmp(nwrdf+i)=zlist(i) 
   enddo
   write(412, rec=recstart) tmp
endif
close(412)
deallocate(tmp)

return
end

subroutine geteridata(nirrep, ilist, filename, isym, numl, numr, nsize)
character*8 filename
integer*8 nsize
integer ilist, nirrep
integer isym, numl(8), numr(8)
!f2py intent(in) :: filename, ilist, nirrep
!f2py intent(out) :: numl(8), numr(8), isym
!f2py intent(out) :: nsize

character*10 filename2
integer*8 irrepx, popp(8), popq(8), popr(8), pops(8), ioppq(8), ioprs(8) 
integer*8 irrep0(16), irrep, rectmp(10,20)
integer*8 moionuml(10,20), moionumr(10,20)
filename2(1:8) = filename
filename2(9:10) =  '.c'

OPEN(UNIT=413, FILE=FILENAME2, FORM='UNFORMATTED')
rewind 413
read(413) IRREPX
read(413) (popp(i), i=1,8)
read(413) (popq(i), i=1,8)
read(413) (popr(i), i=1,8)
read(413) (pops(i), i=1,8)
read(413) (ioppq(i), i = 1,8)
read(413) (ioprs(i), i = 1,8)
read(413) (irrep0(i), i = 1,8)
read(413) (irrep0(i), i = 9,16)
read(413) ((moionuml(i,j),i=1,10),j=1,20)
read(413) ((moionumr(i,j),i=1,10),j=1,20)
read(413) ((rectmp(i,j),i=1,10),j=1,20)
close(413)

isym = irrep0(ilist)

nsize = 0 
do irrep = 1, nirrep 
   nsize = nsize + moionuml(irrep, ilist)*moionumr(irrep, ilist)
enddo
do irrep = 1, 8 
   numl(irrep) = moionuml(irrep, ilist)
   numr(irrep) = moionumr(irrep, ilist) 
enddo
return
end

subroutine getlstso(nirrep, nfirst, nblock, leth, ispin, &
                    ilist, filename, zlist)
implicit none 
integer nfirst, nblock, ispin, ilist, leth
integer nirrep
character*8 filename
real*8  zlist(nblock*leth)

!f2py intent(in) :: nfirst, nblock, ispin, ilist
!f2py intent(in) :: nirrep
!f2py intent(in) :: filename
!f2py intent(out) :: zlist(nblock*leth)
integer*8 irrepx 
integer i, j, k
integer*8 popp(8), popq(8), popr(8), pops(8)
integer*8 ioppq(8), ioprs(8), irrep0(16)
integer*8 moionuml(10,20), moionumr(10,20), rectmp(10,20)
character*10 filename2
real*8, allocatable :: tmp(:)
integer iprcln, iprcwd, nwrdf, nwrdf2
integer*8 nwords, nstart, nend, recstart, recend, irecord 
integer*8 istart, iout, listrec

iprcwd = 4096
iprcln = 4096*8 

filename2(1:8) = filename
filename2(9:10) =  '.c'
OPEN(UNIT=413, FILE=FILENAME2, FORM='UNFORMATTED')
rewind 413
read(413) IRREPX
read(413) (popp(i), i=1,8)
read(413) (popq(i), i=1,8)
read(413) (popr(i), i=1,8)
read(413) (pops(i), i=1,8)
read(413) (ioppq(i), i = 1,8)
read(413) (ioprs(i), i = 1,8)
read(413) (irrep0(i), i = 1,8)
read(413) (irrep0(i), i = 9,16)
read(413) ((moionuml(i,j),i=1,10),j=1,20)
read(413) ((moionumr(i,j),i=1,10),j=1,20)
read(413) ((rectmp(i,j),i=1,10),j=1,20)
close(413)
OPEN(UNIT=412, FILE=FILENAME, FORM='UNFORMATTED', access='direct', recl=iprcln)

listrec = rectmp(ispin, ilist)
nwords = nblock*moionuml(ispin, ilist)
if(nwords<=0) then 
   return
endif
nstart = (nfirst-1)*moionuml(ispin, ilist)
nend   = (nfirst+nblock-1)*moionuml(ispin, ilist)

recstart = (nstart-1)/iprcwd
nwrdf = nstart - iprcwd*recstart 
recstart = recstart + listrec + 1

recend = (nend - 1)/iprcwd
nwrdf2 = nend - iprcwd*recend 
recend = recend + listrec + 1

allocate(tmp(iprcwd))
if(recend>recstart) then 
   read(412, rec=recstart) tmp
   do i = 1, iprcwd-nwrdf  
      zlist(i) = tmp(nwrdf+i)
   enddo
   do iout = 1, recend - recstart - 1 
      irecord = recstart + iout 
      istart = (iout-1)*iprcwd + iprcwd - nwrdf + 1
      read(412, rec=irecord) zlist(istart:istart+iprcwd-1)
   enddo
   read(412, rec=recend) tmp
   istart = nwords - nwrdf2
   do i = 1, nwrdf2 
      zlist(istart + i) = tmp(i) 
   enddo
elseif(recend==recstart) then 
   read(412, rec=recstart) tmp
   do i = 1, nwrdf2 - nwrdf
      zlist(i) = tmp(nwrdf+i)
   enddo
endif
close(412)
deallocate(tmp)

return
end

subroutine save_fock(fijaar, fijaai, fijabr, fijabi, & 
                     fabaar, fabaai, fababr, fababi, &
                     faiaar, faiaai, faiabr, faiabi, & 
                      nva, noa) 
!                     nmo, nva, noa, pop, vrt, isympop, isymvrt) 
!integer nmo, nva, noa, pop(8), vrt(8), isympop(nmo,8), isymvrt(nmo,8)
integer nva, noa
real*8 fijaar(noa,noa), fijaai(noa,noa), fijabr(noa,noa), fijabi(noa,noa) 
real*8 fabaar(nva,nva), fabaai(nva,nva), fababr(nva,nva), fababi(nva,nva) 
real*8 faiaar(nva,noa), faiaai(nva,noa), faiabr(nva,noa), faiabi(nva,noa)

!f2py, intent(in) ::  nva, noa
!f2py, intent(in) ::  fijaar(noa,noa), fijaai(noa,noa), fijabr(noa,noa), fijabi(noa,noa) 
!f2py, intent(in) ::  fabaar(nva,nva), fabaai(nva,nva), fababr(nva,nva), fababi(nva,nva) 
!f2py, intent(in) ::  faiaar(nva,noa), faiaai(nva,noa), faiabr(nva,noa), faiabi(nva,noa)

integer i, j 

OPEN(UNIT=1002, FILE='bdffock.data', FORM='UNFORMATTED')
REWIND 1002
WRITE(1002) ((FIJAAR(I,J), I = 1, NOA), J = 1, NOA)
WRITE(1002) ((FIJAAI(I,J), I = 1, NOA), J = 1, NOA)
WRITE(1002) ((FIJABR(I,J), I = 1, NOA), J = 1, NOA)
WRITE(1002) ((FIJABI(I,J), I = 1, NOA), J = 1, NOA)
WRITE(1002) ((FABAAR(I,J), I = 1, NVA), J = 1, NVA)
WRITE(1002) ((FABAAI(I,J), I = 1, NVA), J = 1, NVA)
WRITE(1002) ((FABABR(I,J), I = 1, NVA), J = 1, NVA)
WRITE(1002) ((FABABI(I,J), I = 1, NVA), J = 1, NVA)
WRITE(1002) ((FAIAAR(I,J), I = 1, NVA), J = 1, NOA)
WRITE(1002) ((FAIAAI(I,J), I = 1, NVA), J = 1, NOA)
WRITE(1002) ((FAIABR(I,J), I = 1, NVA), J = 1, NOA)
WRITE(1002) ((FAIABI(I,J), I = 1, NVA), J = 1, NOA)
CLOSE(1002)
return 
end

subroutine save_symm_uhf(dirprdin, nirrep, nmo, vrt, pop, isymmoin, &
                     symrx, symry, symrz)

integer nmo, symrx, symry, symrz
integer isymmoin(:,:,:)
integer nirrep, pop(8,2), vrt(8,2)
integer dirprdin(8,8)

!f2py, intent(in) :: nmo, symrx, symry, symrz
!f2py, intent(in) :: isymmoin(nmo,8,2)
!f2py, intent(in) :: nirrep, pop(8), vrt(8)
!f2py, intent(in) :: dirprdin(8,8)

integer dirprd(8,8), isymmo(nmo, 8,2)
integer i, j, symmo(8)
integer sx, sy, sz

sx = symrx + 1
sy = symry + 1
sz = symrz + 1

do i = 1, 8
   do j = 1, 8
    dirprd(j,i) = dirprdin(j,i) + 1 
   enddo 
enddo     
         
do i = 1, 8
   symmo(i) = vrt(i,1) + pop(i,1)
enddo        

isymmo=0
do k = 1, 2
do i = 1, nirrep
   do j = 1, pop(i,k)+vrt(i,k)
      isymmo(j,i,k) = isymmoin(j,i,k) + 1
   enddo
enddo
enddo
!write(6,*)'test here'
open(22, file='SYMM_UHF', form='formatted')
rewind 22
write(22,*) nmo
write(22,*) nirrep
write(22,*) dirprd
write(22,*) vrt
write(22,*) pop
write(22,*) symmo
write(22,*) isymmo
write(22,*) sx, sy, sz 
close(22)

return
end

subroutine save_symm(dirprdin, nirrep, nmo, vrt, pop, isymmoin, &
                     symrx, symry, symrz)

integer nmo, symrx, symry, symrz
integer isymmoin(:,:)
integer nirrep, pop(8), vrt(8)
integer dirprdin(8,8)

!f2py, intent(in) :: nmo, symrx, symry, symrz
!f2py, intent(in) :: isymmoin(nmo,8)
!f2py, intent(in) :: nirrep, pop(8), vrt(8)
!f2py, intent(in) :: dirprdin(8,8)

integer dirprd(8,8), isymmo(nmo, 8) 
integer i, j, symmo(8), sx, sy, sz

do i = 1, 8
   do j = 1, 8
    dirprd(j,i) = dirprdin(j,i) + 1
   enddo
enddo

do i = 1, 8
   symmo(i) = vrt(i) + pop(i)
enddo

isymmo=0
do i = 1, nirrep
   do j = 1, pop(i)+vrt(i)
      isymmo(j,i) = isymmoin(j,i) + 1
   enddo
enddo
sx = symrx + 1
sy = symry + 1
sz = symrz + 1
open(22, file='SYMM', form='formatted')
rewind 22
write(22,*) nmo
write(22,*) nirrep
write(22,*) dirprd
write(22,*) vrt
write(22,*) pop
write(22,*) symmo
write(22,*) isymmo
write(22,*) sx, sy, sz 
close(22)

return
end

subroutine mo_pqrs2(iin, dirprdin, nirrep, nmo, numpq, numrs,  &
                   ipospq, iposrs, ioffpqrs,  &
                   popp, popq, popr, pops, &
                   np, nq, nr, ns, &
                   isympin, isymqin, isymrin, isymsin, &
                   ioppr,iopqs, npr, nqs, & 
                   npqrs, filename,ipqrs)
! this is for closed shell situation with eriincore
! to resort (pr|qs) to <pq|rs>
implicit none

character*8 filename
integer nmo, np, nq, nr, ns, ioppr, iopqs, npr, nqs
integer isympin(:,:), isymqin(:,:), isymrin(:,:), isymsin(:,:)
integer nirrep, popp(8), popq(8), popr(8), pops(8)
integer dirprdin(8,8)
integer ipospq(8,8), iposrs(8,8)
integer*8 numpq(8), numrs(8)
integer*8 ioffpqrs(8), npqrs
real*8 iin(:,:)
real*8 ipqrs(npqrs)

!f2py intent(in):: filename
!f2py intent(in):: nmo, np, nq, nr, ns, ioppr, iopqs, npr, nqs
!f2py intent(in):: isympin(nmo,8), isymqin(nmo,8), isymrin(nmo,8), isymsin(nmo,8)
!f2py intent(in):: nirrep, popp(8), popq(8), popr(8), pops(8)
!f2py intent(in):: dirprdin(8,8), numpq(8), numrs(8)
!f2py intent(in):: ipospq(8,8), iposrs(8,8), ioffpqrs(8), npqrs
!f2py intent(in):: iin(npr,nqs)
!f2py intent(out):: ipqrs(npqrs)

integer irreps, irrepq, irrepp, irrepr
integer irreppr, irrepqs, irrepps, irrepqr
integer irreprs, irreppq, irrep, nump, numq, numr, nums
integer irrep1, irrep2, irrep3, irrep4, irrep12, irrep34
integer idxp_in, idxq_in, idxr_in, idxs_in
integer*8 ioffin, ioffout, ioffold
integer p, q, r, s, iinpq, iinrs, itmp, iinpr, iinqs
integer iinqr, iinps, i, j
integer*8 i8
integer dirprd(8,8), isymp(nmo,8), isymq(nmo,8), isymr(nmo,8), isyms(nmo,8)

!write(6,*)'n', npr, np*nr
do i = 1, 8
   do j = 1, 8
    dirprd(j,i) = dirprdin(j,i) + 1
   enddo
enddo


isymp=0; isymq=0; isymr=0; isyms=0
do i = 1, nirrep
   do j = 1, popp(i)
      isymp(j,i) = isympin(j,i) + 1
   enddo
   do j = 1, popq(i)
      isymq(j,i) = isymqin(j,i) + 1
   enddo
   do j = 1, popr(i)
      isymr(j,i) = isymrin(j,i) + 1
   enddo
   do j = 1, pops(i)
      isyms(j,i) = isymsin(j,i) + 1
   enddo
enddo
do j = 1, 8 
! write(6,'(A, 8I4)')'dirprd',  dirprd(j,:)
enddo
do j = 1, nmo
! write(6,'(A, 8I4)')'isymvrta',  isymp(j,:)
enddo
do j = 1, nmo
! write(6,'(A, 8I4)')'isympopa',  isyms(j,:)
enddo
!write(6,*)'pop and vrt', pop, vrt
do irreprs = 1, nirrep
   irreppq = dirprd(irreprs, 1)
   do irreps = 1, nirrep
      irrepr = dirprd(irreps, irreprs)
      numr = popr(irrepr)
      nums = pops(irreps)
      do irrepq = 1, nirrep
         irrepp = dirprd(irrepq, irreppq)
         nump = popp(irrepp)
         numq = popq(irrepq)
         if(ioppr==1.and.iopqs==1) then 
!$omp parallel do
         do q = 1, numq
         do p = 1, nump
         do s = 1, nums
         do r = 1, numr
           idxs_in = isyms(s, irreps)
           idxr_in = isymr(r, irrepr)
           idxq_in = isymq(q, irrepq)
           idxp_in = isymp(p, irrepp)

           iinqs = (idxs_in-1)*nq + idxq_in
           iinpr = (idxr_in-1)*np + idxp_in
           ioffold = ioffout
           ioffout = ioffpqrs(irreprs) - 1 + &
                       (iposrs(irreps, irreprs)-1+(s-1)*numr+r-1)*numpq(irreppq) &
                      + ipospq(irrepq, irreppq) - 1 + (q-1)*nump + p
!          if(ioffout-ioffold/=1) then 
!             write(6,*)'ioffout', ioffout
!          endif
           ipqrs(ioffout) = iin(iinpr, iinqs)
!          if(dabs(ipqrs(ioffout))<1.d-10) then 
!              write(6,'(A,8I9)')'itest', idxp_in, idxq_in, idxr_in, idxs_in, iinpr, iinqs, irreppq, irreprs
!          endif  
         enddo
         enddo
         enddo
         enddo
!$omp end parallel do
         elseif(ioppr==0.and.iopqs==0) then 
!$omp parallel do
         do s = 1, nums
         do r = 1, numr
         do q = 1, numq
         do p = 1, nump
           idxs_in = isyms(s, irreps)
           idxr_in = isymr(r, irrepr)
           idxq_in = isymq(q, irrepq)
           idxp_in = isymp(p, irrepp)

           if(idxp_in>idxr_in) then
              itmp = idxr_in
              idxr_in = idxp_in
              idxp_in = itmp
           endif

           if(idxq_in>idxs_in) then
              itmp = idxs_in
              idxs_in = idxq_in
              idxq_in = itmp
           endif

           iinqs = (idxs_in-1)*idxs_in/2 + idxq_in
           iinpr = (idxr_in-1)*idxr_in/2 + idxp_in
           ioffout = ioffpqrs(irreprs) - 1 + &
                  (iposrs(irreps, irreprs)-1+(s-1)*numr+r-1)*numpq(irreppq)&
                 + ipospq(irrepq, irreppq) - 1 + (q-1)*nump + p
           ipqrs(ioffout) = iin(iinpr, iinqs)
         enddo
         enddo
         enddo
         enddo
!$omp end parallel do
         elseif(ioppr==0.and.iopqs==1) then
!$omp parallel do
         do s = 1, nums
         do r = 1, numr
         do q = 1, numq
         do p = 1, nump
           idxs_in = isyms(s, irreps)
           idxr_in = isymr(r, irrepr)
           idxq_in = isymq(q, irrepq)
           idxp_in = isymp(p, irrepp)

           if(idxp_in>idxr_in) then
              itmp = idxr_in
              idxr_in = idxp_in
              idxp_in = itmp
           endif
           iinqs = (idxs_in-1)*nq + idxq_in
           iinpr = (idxr_in-1)*idxr_in/2 + idxp_in
           ioffout = ioffpqrs(irreprs) - 1 + &
                  (iposrs(irreps, irreprs)-1+(s-1)*numr+r-1)*numpq(irreppq)&
                 + ipospq(irrepq, irreppq) - 1 + (q-1)*nump + p
           ipqrs(ioffout) = iin(iinpr, iinqs)
         enddo
         enddo
         enddo
         enddo
!$omp end parallel do
         elseif(ioppr==1.and.iopqs==0) then
!$omp parallel do
         do s = 1, nums
         do r = 1, numr
         do q = 1, numq
         do p = 1, nump
           idxs_in = isyms(s, irreps)
           idxr_in = isymr(r, irrepr)
           idxq_in = isymq(q, irrepq)
           idxp_in = isymp(p, irrepp)

           if(idxq_in>idxs_in) then
              itmp = idxs_in
              idxs_in = idxq_in
              idxq_in = itmp
           endif

           iinqs = (idxs_in-1)*idxs_in/2 + idxq_in
           iinpr = (idxr_in-1)*np + idxp_in
           ioffout = ioffpqrs(irreprs) - 1 + &
                  (iposrs(irreps, irreprs)-1+(s-1)*numr+r-1)*numpq(irreppq)&
                 + ipospq(irrepq, irreppq) - 1 + (q-1)*nump + p
           ipqrs(ioffout) = iin(iinpr, iinqs)
         enddo
         enddo
         enddo
         enddo
!$omp end parallel do
        endif
       enddo
    enddo
 enddo
!                 write(6,*)'ioffout', ioffout

open(22, file=filename, form='unformatted')
write(22) (ipqrs(i8), i8 = 1, npqrs)
close(22)

return
end 

subroutine getpqrs_incore(dirprdin, nirrep, popp, popq, ioppq, popr, pops, ioprs,  &
                     ioffpqrs, npqrs, ipospq, numpq, iposrs, numrs)
! this is for closed shell situation
! to get the size of eri with p<=q and r<=s 
implicit none
integer dirprdin(8,8)
integer popp(8), popq(8), popr(8), pops(8)
integer ioppq, ioprs
integer ipospq(8,8), iposrs(8,8)
integer nirrep
integer*8 numrs(8), numpq(8)
integer*8 ioffpqrs(8)
integer*8 npqrs
!f2py  intent(in) :: nirrep, ioppq, ioprs
!f2py  intent(in) :: dirprdin(8,8)
!f2py  intent(in) :: popp(8), popq(8), popr(8), pops(8)
!f2py  intent(out):: ioffpqrs(8), ipospq(8,8), iposrs(8,8)
!f2py  intent(out):: numpq(8), numrs(8), npqrs

integer irreps, irrepq, irrepp, irrepr
integer irreprs, irreppq, irrep, nump, numq, numr, nums
integer irrep1, irrep2, irrep3, irrep4, irrep12, irrep34
integer idxp_in, idxq_in, idxr_in, idxs_in
integer nsize, ioffv(8), ioffin, ioffout
integer i, iinpq, iinrs
integer*8 numpq8, numrs8
integer j,  dirprd(8,8)

!write(6,*)'begin getmosize_rhf', nirrep
do i = 1, 8
   do j = 1, 8
    dirprd(j,i) = dirprdin(j,i) + 1
   enddo
enddo
do i = 1, 8
!   write(6,'(A, 8I4)') 'dirprd', dirprd(:,i)
enddo

numpq(1:8) = 0
do irreppq = 1, nirrep
   do irrepq = 1, nirrep
      irrepp = dirprd(irrepq, irreppq)
      nump = popp(irrepp)
      numq = popq(irrepq)
      if(ioppq==1) then 
        numpq(irreppq) = numpq(irreppq) + nump*numq
      elseif(ioppq==0) then 
        if(irrepp>irrepq) cycle
        if(irrepp==irrepq) then
           numpq(irreppq) = numpq(irreppq) + nump*(nump-1)/2
        elseif(irrep<irrepq) then
           numpq(irreppq) = numpq(irreppq) + nump*numq
        endif
      elseif(ioppq==2) then 
        if(irrepp>irrepq) cycle
        if(irrepp==irrepq) then
           numpq(irreppq) = numpq(irreppq) + nump*(nump+1)/2
        elseif(irrep<irrepq) then
           numpq(irreppq) = numpq(irreppq) + nump*numq
        endif
      endif
   enddo
enddo

numrs(1:8) = 0
do irreppq = 1, nirrep
   do irrepq = 1, nirrep
      irrepp = dirprd(irrepq, irreppq)
      nump = popr(irrepp)
      numq = pops(irrepq)
      if(ioprs==1) then
        numrs(irreppq) = numrs(irreppq) + nump*numq
      elseif(ioprs==0) then
        if(irrepp>irrepq) cycle
        if(irrepp==irrepq) then
           numrs(irreppq) = numrs(irreppq) + nump*(nump-1)/2
        elseif(irrep<irrepq) then
           numrs(irreppq) = numrs(irreppq) + nump*numq
        endif
      elseif(ioprs==2) then
        if(irrepp>irrepq) cycle
        if(irrepp==irrepq) then
           numrs(irreppq) = numrs(irreppq) + nump*(nump+1)/2
        elseif(irrep<irrepq) then
           numrs(irreppq) = numrs(irreppq) + nump*numq
        endif
      endif
   enddo
enddo

do irreprs = 1, nirrep
   iposrs(1,irreprs) = 1
   do irreps = 1, nirrep-1
      irrepr = dirprd(irreprs, irreps) 
      nums = pops(irreps) 
      numr = popr(irrepr)    
      if(ioprs==1) then 
        iposrs(irreps+1, irreprs) = iposrs(irreps, irreprs) + nums*numr
      elseif(ioprs==0) then 
        if(irrepr>irreps) then 
           iposrs(irreps+1, irreprs) = iposrs(irreps, irreprs)
        elseif(irrepr==irreps) then 
           iposrs(irreps+1, irreprs) = iposrs(irreps, irreprs) + nums*(nums-1)/2
        elseif(irrepr<irreps) then 
           iposrs(irreps+1, irreprs) = iposrs(irreps, irreprs) + nums*numr
        endif
      elseif(ioprs==2) then 
        if(irrepr>irreps) then 
           iposrs(irreps+1, irreprs) = iposrs(irreps, irreprs)
        elseif(irrepr==irreps) then 
           iposrs(irreps+1, irreprs) = iposrs(irreps, irreprs) + nums*(nums+1)/2
        elseif(irrepr<irreps) then 
           iposrs(irreps+1, irreprs) = iposrs(irreps, irreprs) + nums*numr
        endif
      endif
   enddo
enddo

do irreprs = 1, nirrep
   ipospq(1,irreprs) = 1
   do irreps = 1, nirrep-1
      irrepr = dirprd(irreprs, irreps)
      nums = popq(irreps)
      numr = popp(irrepr)
      if(ioprs==1) then
        ipospq(irreps+1, irreprs) = ipospq(irreps, irreprs) + nums*numr
      elseif(ioprs==0) then
        if(irrepr>irreps) then
           ipospq(irreps+1, irreprs) = ipospq(irreps, irreprs)
        elseif(irrepr==irreps) then
           ipospq(irreps+1, irreprs) = ipospq(irreps, irreprs) + nums*(nums-1)/2
        elseif(irrepr<irreps) then
           ipospq(irreps+1, irreprs) = ipospq(irreps, irreprs) + nums*numr
        endif
      elseif(ioprs==2) then
        if(irrepr>irreps) then
           ipospq(irreps+1, irreprs) = ipospq(irreps, irreprs)
        elseif(irrepr==irreps) then
           ipospq(irreps+1, irreprs) = ipospq(irreps, irreprs) + nums*(nums+1)/2
        elseif(irrepr<irreps) then
           ipospq(irreps+1, irreprs) = ipospq(irreps, irreprs) + nums*numr
        endif
      endif
   enddo
enddo

ioffpqrs(1) = 1
do irreprs = 1, nirrep-1
   irreppq = dirprd(irreprs, 1)
   ioffpqrs(irreprs+1) = ioffpqrs(irreprs) + numpq(irreppq)*numrs(irreprs)
enddo

npqrs = 0
do irreprs = 1, nirrep
   irreppq = dirprd(irreprs, 1)
   npqrs = npqrs + numpq(irreppq)*numrs(irreprs)
enddo
!write(6,'(A,8I5)')'symmo', symmo(:)
!write(6,'(A,8I5)')'symmo2', symmo2(:)
!write(6,'(A,8I8)')'ioffmo', ioffmo(:)
!write(6,*)'size of pqrs with symmetry', mosize

 return
 end

subroutine idsymsoc(dirprdin, nirrep, popp, popq, ioppq, popr, pops, ioprs,  &
                      npqrs)
! this is for closed shell situation
! to get the size of eri with p<=q and r<=s 
implicit none
integer dirprdin(8,8)
integer popp(8), popq(8), popr(8), pops(8)
integer ioppq, ioprs
integer nirrep
integer*8 npqrs
!f2py  intent(in) :: nirrep, ioppq, ioprs
!f2py  intent(in) :: dirprdin(8,8)
!f2py  intent(in) :: popp(8), popq(8), popr(8), pops(8)
!f2py  intent(out):: npqrs

integer irreps, irrepq, irrepp, irrepr
integer irreprs, irreppq, irrep, nump, numq, numr, nums
integer irrep1, irrep2, irrep3, irrep4, irrep12, irrep34
integer idxp_in, idxq_in, idxr_in, idxs_in
integer nsize, ioffv(8), ioffin, ioffout
integer i, iinpq, iinrs
integer*8 numpq(8), numrs(8)
integer j,  dirprd(8,8)

!write(6,*)'begin getmosize_rhf', nirrep
do i = 1, 8
   do j = 1, 8
    dirprd(j,i) = dirprdin(j,i) + 1
   enddo
enddo
do i = 1, 8
!   write(6,'(A, 8I4)') 'dirprd', dirprd(:,i)
enddo

numpq(1:8) = 0
do irreppq = 1, nirrep
   do irrepq = 1, nirrep
      irrepp = dirprd(irrepq, irreppq)
      nump = popp(irrepp)
      numq = popq(irrepq)
      if(ioppq==1) then 
        numpq(irreppq) = numpq(irreppq) + nump*numq
      elseif(ioppq==0) then 
        if(irrepp>irrepq) cycle
        if(irrepp==irrepq) then
           numpq(irreppq) = numpq(irreppq) + nump*(nump-1)/2
        elseif(irrep<irrepq) then
           numpq(irreppq) = numpq(irreppq) + nump*numq
        endif
      elseif(ioppq==2) then 
        if(irrepp>irrepq) cycle
        if(irrepp==irrepq) then
           numpq(irreppq) = numpq(irreppq) + nump*(nump+1)/2
        elseif(irrep<irrepq) then
           numpq(irreppq) = numpq(irreppq) + nump*numq
        endif
      endif
   enddo
enddo

numrs(1:8) = 0
do irreppq = 1, nirrep
   do irrepq = 1, nirrep
      irrepp = dirprd(irrepq, irreppq)
      nump = popr(irrepp)
      numq = pops(irrepq)
      if(ioprs==1) then
        numrs(irreppq) = numrs(irreppq) + nump*numq
      elseif(ioprs==0) then
        if(irrepp>irrepq) cycle
        if(irrepp==irrepq) then
           numrs(irreppq) = numrs(irreppq) + nump*(nump-1)/2
        elseif(irrep<irrepq) then
           numrs(irreppq) = numrs(irreppq) + nump*numq
        endif
      elseif(ioprs==2) then
        if(irrepp>irrepq) cycle
        if(irrepp==irrepq) then
           numrs(irreppq) = numrs(irreppq) + nump*(nump+1)/2
        elseif(irrep<irrepq) then
           numrs(irreppq) = numrs(irreppq) + nump*numq
        endif
      endif
   enddo
enddo

npqrs = 0
do irreprs = 1, nirrep
   irreppq = dirprd(irreprs, 1)
   npqrs = npqrs + numpq(irreppq)*numrs(irreprs)
!   npqrs = npqrs + numpq(irreppq)*numrs(irreprs)
enddo

 return
 end

subroutine mo_pqrs_sym(iin, dirprdin, nirrep, nmo, numpq, numrs,  &
                   ipospq, iposrs, ioffpqrs,  &
                   popp, popq, popr, pops, &
                   np, nq, nr, ns, &
                   isympin, isymqin, isymrin, isymsin, &
                   ioppr,iopqs, npr, nqs, &
                   npqrs, filename,ipqrs)
! this is for closed shell situation with eriincore
! to resort <pq|rs>+<pq|sr> = (pr|qs)+(ps|qr)
implicit none

character*8 filename
integer nmo, np, nq, nr, ns, ioppr, iopqs, npr, nqs
integer isympin(:,:), isymqin(:,:), isymrin(:,:), isymsin(:,:)
integer nirrep, popp(8), popq(8), popr(8), pops(8)
integer dirprdin(8,8)
integer ipospq(8,8), iposrs(8,8)
integer*8 numpq(8), numrs(8)
integer*8 ioffpqrs(8), npqrs
real*8 iin(:,:)
real*8 ipqrs(npqrs)

!f2py intent(in):: filename
!f2py intent(in):: nmo, np, nq, nr, ns, ioppr, iopqs, npr, nqs
!f2py intent(in):: isympin(nmo,8), isymqin(nmo,8), isymrin(nmo,8),
!isymsin(nmo,8)
!f2py intent(in):: nirrep, popp(8), popq(8), popr(8), pops(8)
!f2py intent(in):: dirprdin(8,8), numpq(8), numrs(8)
!f2py intent(in):: ipospq(8,8), iposrs(8,8), ioffpqrs(8), npqrs
!f2py intent(in):: iin(npr,nqs)
!f2py intent(out):: ipqrs(npqrs)

integer irreps, irrepq, irrepp, irrepr
integer irreppr, irrepqs, irrepps, irrepqr
integer irreprs, irreppq, irrep, nump, numq, numr, nums
integer irrep1, irrep2, irrep3, irrep4, irrep12, irrep34
integer idxp_in, idxq_in, idxr_in, idxs_in
integer*8 ioffin, ioffout, ioffold
integer p, q, r, s, iinpq, iinrs, itmp, iinpr, iinqs
integer iinqr, iinps, i, j
integer*8 i8
integer dirprd(8,8), isymp(nmo,8), isymq(nmo,8), isymr(nmo,8), isyms(nmo,8)

!write(6,*)'n', npr, np*nr
do i = 1, 8
   do j = 1, 8
    dirprd(j,i) = dirprdin(j,i) + 1
   enddo
enddo


isymp=0; isymq=0; isymr=0; isyms=0
do i = 1, nirrep
   do j = 1, popp(i)
      isymp(j,i) = isympin(j,i) + 1
   enddo
   do j = 1, popq(i)
      isymq(j,i) = isymqin(j,i) + 1
   enddo
   do j = 1, popr(i)
      isymr(j,i) = isymrin(j,i) + 1
   enddo
   do j = 1, pops(i)
      isyms(j,i) = isymsin(j,i) + 1
   enddo
enddo
do j = 1, 8
! write(6,'(A, 8I4)')'dirprd',  dirprd(j,:)
enddo
do j = 1, nmo
! write(6,'(A, 8I4)')'isymvrta',  isymp(j,:)
enddo
do j = 1, nmo
! write(6,'(A, 8I4)')'isympopa',  isyms(j,:)
enddo
!write(6,*)'pop and vrt', pop, vrt
do irreprs = 1, nirrep
   irreppq = dirprd(irreprs, 1)
   do irreps = 1, nirrep
      irrepr = dirprd(irreps, irreprs)
      if(irrepr>irreps) cycle
      numr = popr(irrepr)
      nums = pops(irreps)
      do irrepq = 1, nirrep
         irrepp = dirprd(irrepq, irreppq)
         if(irrepp>irrepq) cycle
         nump = popp(irrepp)
         numq = popq(irrepq)
         if(irrepr<irreps.and.irrepp<irrepq) then
!$omp parallel do
         do q = 1, numq
         do p = 1, nump
         do s = 1, nums
         do r = 1, numr
           idxs_in = isyms(s, irreps)
           idxr_in = isymr(r, irrepr)
           idxq_in = isymq(q, irrepq)
           idxp_in = isymp(p, irrepp)

           iinqs = (idxs_in-1)*nq + idxq_in
           iinpr = (idxr_in-1)*np + idxp_in
           ioffold = ioffout
           ioffout = ioffpqrs(irreprs) - 1 + &
                       (iposrs(irreps, irreprs)-1+(s-1)*numr+r-1)*numpq(irreppq)&
                      + ipospq(irrepq, irreppq) - 1 + (q-1)*nump + p
!          if(ioffout-ioffold/=1) then 
!             write(6,*)'ioffout', ioffout
!          endif
           ipqrs(ioffout) = iin(iinpr, iinqs)
!          if(dabs(ipqrs(ioffout))<1.d-10) then 
!              write(6,'(A,8I9)')'itest', idxp_in, idxq_in, idxr_in, idxs_in,
!              iinpr, iinqs, irreppq, irreprs
!          endif  
         enddo
         enddo
         enddo
         enddo
!$omp end parallel do
         elseif(ioppr==0.and.iopqs==0) then
!$omp parallel do
         do s = 1, nums
         do r = 1, numr
         do q = 1, numq
         do p = 1, nump
           idxs_in = isyms(s, irreps)
           idxr_in = isymr(r, irrepr)
           idxq_in = isymq(q, irrepq)
           idxp_in = isymp(p, irrepp)

           if(idxp_in>idxr_in) then
              itmp = idxr_in
              idxr_in = idxp_in
              idxp_in = itmp
           endif

           if(idxq_in>idxs_in) then
              itmp = idxs_in
              idxs_in = idxq_in
              idxq_in = itmp
           endif

           iinqs = (idxs_in-1)*idxs_in/2 + idxq_in
           iinpr = (idxr_in-1)*idxr_in/2 + idxp_in
           ioffout = ioffpqrs(irreprs) - 1 + &
                  (iposrs(irreps, irreprs)-1+(s-1)*numr+r-1)*numpq(irreppq)&
                 + ipospq(irrepq, irreppq) - 1 + (q-1)*nump + p
           ipqrs(ioffout) = iin(iinpr, iinqs)
         enddo
         enddo
         enddo
         enddo
!$omp end parallel do
         elseif(ioppr==0.and.iopqs==1) then
!$omp parallel do
         do s = 1, nums
         do r = 1, numr
         do q = 1, numq
         do p = 1, nump
           idxs_in = isyms(s, irreps)
           idxr_in = isymr(r, irrepr)
           idxq_in = isymq(q, irrepq)
           idxp_in = isymp(p, irrepp)

           if(idxp_in>idxr_in) then
              itmp = idxr_in
              idxr_in = idxp_in
              idxp_in = itmp
           endif
           iinqs = (idxs_in-1)*nq + idxq_in
           iinpr = (idxr_in-1)*idxr_in/2 + idxp_in
           ioffout = ioffpqrs(irreprs) - 1 + &
                  (iposrs(irreps, irreprs)-1+(s-1)*numr+r-1)*numpq(irreppq)&
                 + ipospq(irrepq, irreppq) - 1 + (q-1)*nump + p
           ipqrs(ioffout) = iin(iinpr, iinqs)
         enddo
         enddo
         enddo
         enddo
!$omp end parallel do
         elseif(ioppr==1.and.iopqs==0) then
!$omp parallel do
         do s = 1, nums
         do r = 1, numr
         do q = 1, numq
         do p = 1, nump
           idxs_in = isyms(s, irreps)
           idxr_in = isymr(r, irrepr)
           idxq_in = isymq(q, irrepq)
           idxp_in = isymp(p, irrepp)

           if(idxq_in>idxs_in) then
              itmp = idxs_in
              idxs_in = idxq_in
              idxq_in = itmp
           endif

           iinqs = (idxs_in-1)*idxs_in/2 + idxq_in
           iinpr = (idxr_in-1)*np + idxp_in
           ioffout = ioffpqrs(irreprs) - 1 + &
                  (iposrs(irreps, irreprs)-1+(s-1)*numr+r-1)*numpq(irreppq)&
                 + ipospq(irrepq, irreppq) - 1 + (q-1)*nump + p
           ipqrs(ioffout) = iin(iinpr, iinqs)
         enddo
         enddo
         enddo
         enddo
!$omp end parallel do
        endif
       enddo
    enddo
 enddo
!                 write(6,*)'ioffout', ioffout

open(22, file=filename, form='unformatted')
write(22) (ipqrs(i8), i8 = 1, npqrs)
close(22)

return
end
 
